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Abstract 

Natural selection drives populations towards higher fitness, but crossing fitness valleys or plateaus may facilitate progress 
up a rugged fitness landscape involving epistasis. We investigate quantitatively the effect of subdividing an asexual 
population on the time it takes to cross a fitness valley or plateau. We focus on a generic and minimal model that includes 
only population subdivision into equivalent demes connected by global migration, and does not require significant size 
changes of the demes, environmental heterogeneity or specific geographic structure. We determine the optimal speedup of 
valley or plateau crossing that can be gained by subdivision, if the process is driven by the deme that crosses fastest. We 
show that isolated demes have to be in the sequential fixation regime for subdivision to significantly accelerate crossing. 
Using IVlarkov chain theory, we obtain analytical expressions for the conditions under which optimal speedup is achieved: 
valley or plateau crossing by the subdivided population is then as fast as that of its fastest deme. We verify our analytical 
predictions through stochastic simulations. We demonstrate that subdivision can substantially accelerate the crossing of 
fitness valleys and plateaus in a wide range of parameters extending beyond the optimal window. We study the effect of 
varying the degree of subdivision of a population, and investigate the trade-off between the magnitude of the optimal 
speedup and the width of the parameter range over which it occurs. Our results, obtained for fitness valleys and plateaus, 
also hold for weakly beneficial intermediate mutations. Finally, we extend our work to the case of a population connected 
by migration to one or several smaller islands. Our results demonstrate that subdivision with migration alone can 
significantly accelerate the crossing of fitness valleys and plateaus, and shed light onto the quantitative conditions 
necessary for this to occur. 
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Introduction 

Natural selection drives populations towards higher fitness (i.e. 
reproductive success), but crossing fitness valleys or plateaus may 
facilitate progress up a rugged fitness landscape. Rugged fitness 
landscapes arise from epistasis, i.e. interactions between genetic 
variants. For instance, two mutations together can yield a benefit 
while each of them alone is detrimental: such reciprocal sign 
epistasis can give rise to a fitness valley [1,2]. While the high 
dimensionality of genotype space makes it challenging to probe the 
structure of fitness landscapes [3,4], evidence has been accumu- 
lating for frequent landscape ruggedness, especially in recent years 
[1,2,4-1,5]. 

Population structure can play an important role in evolution 
[16-24]. In particular, the time taken to cross a fitness valley or 
plateau depends on population size since stochastic effects such as 
genetic drift have an increased importance in small populations, 
allowing neutral and deleterious mutations to fix with increased 
probability [25-28]. Population subdivision into demes can allow 
the maintenance of larger genetic diversity due to increased 
genetic drift as well as to the quasi-independent explorations of the 



fitness landscape that are run in parallel by each deme. 
Subdivision may thereby facilitate valley or plateau crossing 
locally and subsequent migration can then spread beneficial 
mutations throughout the entire subdivided population ("meta- 
population"). This idea was first discussed by Wright in his shifting 
balance theory [29-32] and the importance of this effect has been 
the subject of a long debate [33-42]. In this work, we investigate 
the role of subdivision with global migration alone, without 
additional effects such as strong dependence of deme size on 
fitness, including extinction and refounding of demes, which 
played a crucial role in Wright's theory. Our generic and minimal 
model enables us to quantatively determine the conditions under 
which population subdivision accelerates fitness valley or plateau 
crossing. 

Studying quantitatively the effect of subdivision on evolution 
may help in inferring fitness landscape structure from evolution 
experiments [43] . Work on structured populations has been used 
as qualitative proof of landscape ruggedness [16]. Current 
experiments investigating the evolution of subdivided populations 
at various migration rates have produced mixed results, some 
demonstrating faster adaptation of subdivided populations [44], 
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Author Summary 

Experimental evidence has recently been accumulating to 
suggest that fitness landscape ruggedness Is common in a 
variety of organisms. Rugged landscapes arise from 
interactions between genetic variants, called epistasis, 
which can lead to fitness valleys or plateaus. The time 
needed to cross such fitness valleys or plateaus exhibits a 
rich dependence on population size, since stochastic 
effects have higher Importance in small populations, 
increasing the probability of fixation of neutral or 
deleterious mutants. This may lead to an advantage of 
population subdivision, a possibility which has been 
strongly debated for nearly one hundred years. In this 
work, we quantitatively determine when, and to what 
extent, population subdivision accelerates valley and 
plateau crossing. Using the simple model of an asexual 
population subdivided Into Identical demes connected by 
gobal migration, we derive the conditions under which 
crossing by a subdivided population is driven by its fastest 
deme, thus giving rise to the maximal speedup. Our 
analytical predictions are verified using stochastic simula- 
tions. We Investigate the effect of varying the degree of 
subdivision of a population. We generalize our results to 
weakly beneficial Intermediates and to different popula- 
tion structures. We discuss the magnitude and robustness 
of the effect for realistic parameter values. 

and others not [45]. It is therefore important to determine under 
what conditions subdivision accelerates fitness valley or plateau 
crossing. Additionally, population subdivision is extremely com- 
mon in natural systems. For instance, evidence has recently been 
found for compartmentalization of HIV in different organs of a 
single patient [46,47]. 

Here we show that subdivision can sigiiificandy accelerate 
fitness valley or plateau crossing over a wide parameter range, 
both with respect to a non-subdivided population and with respect 
to a single deme. Intuitively, deleterious or neutral intermediate 
mutations may fix in individual demes, allowing for the 
maintenance of a larger proportion of these mutants in a 
metapopulation than in a well-mixed population. We fu-st 
determine the optimal speedup of valley or plateau crossing by 
subdivision, in the best possible scenario where valley or plateau 
crossing by the metapopulation is driven by that of its fastest deme. 
This enables us to demonstrate that isolated demes must be in the 
sequential frxation regime for subdivision to significantly accelerate 
crossing. We then determine the conditions under which the best 
possible scenario can be realized. Using Markov chain theory, we 
obtain analytical expressions for the parameter range where valley 
or plateau crossing by a metapopulation is as fast as that of its 
fastest deme. Our analytical predictions are verified using 
stochastic simulations. Furthermore, we discuss the effect of 
varying the degree of subdivision of a population, and investigate 
the trade-off between the magnitude of the optimal speedup and 
the width of the parameter range over which it occurs. Finally, we 
extend our work to weakly beneficial mutations and to a 
population connected to smaller islands, and we discuss the 
magnitude and robustness of the effect for realistic parameter 
values. 

Results 

Our results are organized as follows. First, we specify our model 
for the evolutionary dynamics of a subdivided population with 
migration. Then, we focus on the 'best possible' scenario where the 



metapopulation is driven by its fastest deme. We calculate the ratio 
of the valley-crossing time for the metapopulation to the valley- 
crossing time for an equally-sized well-mixed population under 
this strong assumption. This yields the optimal speedup that may 
be obtained by subdivision, and enables us to demonstrate that 
sequential fixation in individual demes is necessary to achieve a 
significant speedup. Then, we determine the range of parameter 
values for which the best possible scenario is attained, i.e. the 
valley-crossing time for the metapopulation is indeed dominated 
by the valley-crossing time of its fastest deme. Qualitatively, 
migration has to be both rare enough to enable demes to cross the 
fitness valley or plateau quasi-independently and frequent enough 
to allow fast spreading of the final beneficial mutation to the whole 
metapopulation once it has frxed in the fastest deme: these 
conditions yield an optiinal window of migration rates. Finally, we 
compare our analytical predictions with results from stochastic 
simulations. 

Model of evolutionary dynamics in a subdivided 
population 

We focus on asexual individuals, characterized by their 
genotype and associated fitness /. Each individual has a division 
rate proportional to /, and a death rate d, which is the same for 
all. We consider an ensemble of D identical demes, each with a 
constant number of individuals. The division rate averaged over 
the individuals of a deme is thus equal to the death rate d. We 
treat migration as a random exchange of two individuals between 
two different demes, occurring at rate 2m per individual. In our 
model, exchange between any two demes is equally likely, as in 
Wright's "island model" [29]. This constitutes a generic and 
minimal model of subdivision with migration, without any 
dependence of migration rate on the average fitness of a deme 
(in contrast with models where demes containing beneficial 
mutants increase significantly in size and migrate more rapidly 
[30,33]), or additional effects of extinction and re-founding of 
demes [30,32,33], specific geographic structure [16,17,19-21], or 
spatially heterogeneous environments [18,22-24], on which 
previous studies focused. 

We consider the simplest fitness valley or plateau, involving 
three successive genotypes denoted by '0', '1' and '2' (see Fig. lA). 
The initial genotype is taken as reference for fitness: yo = l- We 
denote the fitnesses of the subsequent genotypes hy f\ = \—& and 
f2 = \-\-s. The first mutation is assumed to be either neutral 
((5 = 0), which yields a fitness plateau, or deleterious [5 > 0), which 
corresponds to a fitness valley, while the second mutation is 
assumed to be beneficial {s>Q). We focus on first mutations that 
are not too strongly deleterious: S«l. We only allow forward 
mutations, and note that including back mutations does not 
qualitatively affect crossing times [28] . Finally, we assume that all 
mutations have probability /J per division, but generalization to 
different mutation probabilities is straightforward. 

In this paper, we focus on the average time T„, required for the 
whole metapopulation to cross the fitness valley or plateau, i.e. to 
fix mutation '2' in all demes, starting from an initial state where all 
individuals have genotype '0'. 

The best possible scenario 

For small enough migration rates, each deme in the metapop- 
ulation performs a quasi-independent trial at crossing the valley or 
plateau. At best, the valley or plateau crossing time T„, of the whole 
metapopulation is dominated by that, T^, of the "champion" deme 
in the metapopulation, i.e. the deme that crosses the fitness valley 
or plateau fastest. 
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Figure 1. Population subdivision with migration can accelerate fitness valley crossing. A. Fitness valley: fitness / versus genotype g. B. 
Schematic representation of the best possible scenario, for a metapopulation with £) = 7 demes. Each square represents a deme of identical size, and 
a row represents the metapopulation. Colors represent genotypes, with the color-code defined in A. Initially (top row), all demes have genotype '0'. 
The demes explore the fitness landscape described in A quasi-independently, and one of them, the "champion" deme (second from the left here), 
crosses the fitness valley first (second and third row). Individual demes are assumed to be in the sequential fixation regime, so this deme fixes first 
mutation '1' and then mutation '2'. The beneficial mutation '2' then spreads by migration, which is modeled by random exchange of individuals 
between demes (arrow on the fourth row), leading to fixation of mutation '2' in the whole metapopulation (fifth row). C. Average valley crossing time 
T of a non-structured population, as a function of its carrying capacity K, in logarithmic scale. Dots are simulation results, averaged over 1 000 runs for 
each value of K; error bars represent 95% confidence intervals (CI). Theoretical predictions from Ref. [28] are plotted for the sequential fixation regime 
(blue line) and for the tunneling regime (red line), using N = 0.9K (see text) to make the correspondence. The transition between these two regimes 
is indicated by a dotted line. The carrying capacities at stake in D are highlighted in green {id: isolated deme; ns: non-subdivided population). 
Parameter values: (/ = 0.1, /i = 8 x 10^', s = 0.3 and <5 = 6xlO"'. D. Average valley crossing time t„, of a metapopulation composed oi D = l demes 
each with carrying capacity = 357 (total carrying capacity: DK = 2A99), plotted versus the migration-to-mutation rate ratio m/(pid), in logarithmic 
scale. Parameter values are the same as in C, and only the migration rate m is varied. Dots represent simulation results averaged over 1000 runs for 
each value of m, and error bars are 95% CI. Black vertical lines represent the limits of the interval of m/(nd) in Eq. 14. Blue (resp. red) line: valley 
crossing time for an isolated deme (id) with K = 351 (resp. a non-subdivided population (ns) with ^ = 2500) for the same parameter values, averaged 
over 1000 runs; shaded regions: 95% CI. Dashed blue (resp. red) lines: corresponding theoretical predictions from Ref. [28] (see C). 
doi:10.1371/journal.pcbi.1003778.g001 



We now focus on this best possible scenario, wliich is illustrated 
schematically in Fig. IB: first, the champion deme crosses the 
valley or plateau by sequential fixation, and then the beneficial 
mutation rapidly spreads by migration of through the whole 
metapopulation. Once this best possible scenario is characterized, 
the crucial question will be whether, and under what conditions, it 
can be attained: this point will be addressed in the following 
section. 

Determination of T^. Valley or plateau crossing by a non- 
structured, well-mixed population can occur by two different 
mechanisms: sequential fixation and tunneling. The former 
corresponds to fixation of mutation ' 1' in the whole population, 
and to subsequent frxation of the beneficial mutation '2'. 
Conversely, the latter occurs when the beneficial mutation arises 
in a small fluctuating minority of first-mutants, and fnces directly: 
tunneling does not involve fixation of the intermediate mutation 
'1' [28]. For given values of the parameters S, s, and fi, sequential 
fixation is the fastest process for small populations, where genetic 
drift plays an important part. Tunneling becomes the dominant 
process of vaUey or plateau crossing when the number N of 
individuals per deme exceeds a threshold value Nx , which 
depends on 5, s, and ji (see Ref [28] for a full discussion of this 
threshold value). Fig. IC shows simulation results for the valley 
crossing time T of a non-subdivided population versus its size, and 
illustrates these two different regimes and the transition between 
them. Note that in our simulations (described in Methods, Sec. 1), 
we hold fixed the carrying capacity K of populations (or demes) 
instead of the number of individuals A^. This softer constraint is 
more realistic and avoids some possible biases in the metapopu- 
lation case (see Methods, Sec. 1.2). In practice, each individual 



divides at a rate /(I —N/K) and dies at a constant rate d: hence, 
at steady-state, N x K{1 —d /f). We choose d = 0.\, and fitnesses / 
of order one, thus Nx0.9K. 

We now consider D independent demes with no migration, and 
we determine the crossing time of the fastest of these D demes, 
both for demes in the sequential fixation regime and for demes in 
the tunneling regime. 

Demes in the sequential fixation regime. Let 

denote the probability of fixation of genotype 'f, with fitness fj, 
starting from a single individual with genotype 'f in a deme where 
all other individuals initially have genotype 'f and fitness fi¥^fj 
[25,28]. If fi=fj, the probability of fixation of genotype '/ reads 
Pij = l/N. Valley or plateau crossing by sequential fixation 
involves two successive steps. The first step, fixation of the 
intermediate mutation '1', occurs with rate r^i = N f.idpoi , where 
Nfid is the total mutation rate in the deme. (Recall that the deme 
size is fixed, and that d represents the birth/death rate. Note 
that the correspondence with Ref [28] is obtained by multiplying 
by 1 /d all the timescales in this reference, which are expressed in 
numbers of generations.) Similarly, the second step, fixation of the 
final beneficial mutation '2', has rate rn =Njidp\2. The first step is 
longer than the second one since mutation '1' is neutral or 
deleterious, while mutation '2' is beneficial. If the first step 
dominates, the distribution of crossing times is approximately 
exponential with rate roi. The shortest crossing time among D 
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independent demes is then distributed exponentially with rate 
-Droi (see Methods, Sec. 2). Thus, the average crossing time of the 



champion deme reads Tc«(-Dfoi) . Denoting by t,-^? 
average crossing time for an isolated deme, we obtain 



'01 



1 

D 



the 



(2) 



Hence, the champion deme crosses the valley D times faster on 
average than a single deme. This simple result holds for 
Dpm«p\2. For simplicity, we restrict ourselves to this regime in 
the main text, but we provide the general method for calculating 
Tc in Methods, Sec. 2. We use this general method to calculate 
numerically the exact value of t^. in our examples below. 

Demes in the tunneling regime. Assuming that A^/i < 1, so 
that there is no competition between different mutant lineages, 
valley or plateau crossing by tunneling involves a single event with 
constant rate, namely the appearance of a "successful" 'I'-mutant, 
whose lineage includes a '2'-mutant that fixes [28]. Crossing time 
is thus exponentially distributed. Therefore, in this case too, the 
crossing time 1^ of the champion deme among D isolated demes is 
D times smaller than that of an average isolated deme (see 
Methods, Sec. 2): Eq. 2 is valid in the tunneling regime too. 

Sequential fixation in individual demes is necessary for 
significant speedups. In the best possible scenario, where 
crossing by the metapopulation is dominated by that of the 
champion deme, i.e. Tm~Tc, the previous paragraph shows that 
Tm/Taa;l/Z), both when isolated demes are in the sequential 
fixation regime and when they are in the tunneling regime. Hence, 
it is necessary to have 



(3) 



where t„ is the average crossing time of the non-subdivided 
population, for subdivision to speed up valley or plateau crossing 
in the best scenario (i.e. for t„ « to be smaller than t„s). This 
necessary condition is general since it holds a fortiori beyond the 
best scenario. Graphically, in Fig. IC, which is a logarithmic plot 
of crossing time versus population size for a non-structured 
population, the slope of the line joining the isolated deme to the 
non-subdivided population has to be less negative than -1 in order 
for speedups to be possible. Recall indeed that the nonsubdivided 
population is D times larger than an isolated deme. The necessary 
condition in Eq. 3 leaves the possibility of significant speedups in 
the non-trivial case where a single isolated deme crosses slower 
than a non-subdivided population (Tut > t„s). Fig. ID demonstrates 
a significant speedup by subdivision obtained in this regime where 

l<Ta/T„.5<Z). 

Let us [:onsider a metapopulation such that isolated demes are 
in the tunneling regime. Then, the larger non-subdivided 
population with ND individuals is also in the tunneling regime 
[28]. Assuming that NDn<\, valley or plateau crossing by this 
non-subdivided population follows the same laws as crossing by 
the demes. Since the average crossing time by tunneling is 
inversely proportional to population size (see Ref [28] and 
Fig. IC), we obtain Ta/T„,5 = Z), in contradiction with Eq. 3. This 
implies that, even in the best possible scenario, subdivision cannot 
accelerate crossing if isolated demes are in the tunneling regime 
(since here, t,„/t„5«1). Thus, having isolated demes in the 
sequential fixation regime is a necessary condition for subdivision 
to accelerate crossing. Importantly, however, the non-subdivided 



population is not required to be in the sequential fixation regime. 
For instance, in Fig. ID, the non-subdivided population is in the 
tunneling regime. Note that when ND/i> 1, the, population enters 
the semi-deterministic regime [28] and the average crossing time 
need not be proportional to \/N. Minor speedups may exist in this 
regime, but such effects are beyond the scope of this work. In all 
the following, we wiU focus on the regime NDfi< 1. 

Maximal possible speedup by subdivision. The speedup 
gained by subdividing a population of a given total size is directly 
described by the ratio tm/t„s of the valley crossing time of a 
metapopulation to that of a non-subdivided population. Here, we 
discuss the values this ratio can take in the best possible scenario, 
where valley crossing by the metapopulation is dominated by that 
of the champion deme, and we determine the valley depth for 
which the highest speedups are obtained (i.e. for which this ratio is 
smallest). 

Let us first focus on the case where both the non-subdivided 
population and the isolated deme are in the sequential fixation 
regime. The average valley crossing time by the champion deme 
reads ZcXl/(DNfidpQ\) (see our calculation of above). In the 
best possible scenario, T„, xt^.. The average valley crossing time by 
the non-subdivided population is T„s^^/(DNfidp'oi), where 
P'oi =(fi^ — — 1) is the fixation probability of an individ- 

ual with genotype ' 1' in a population of ND individuals where all 
the others initially have genotype '0' (see Eq. 1). Hence, we obtain 



' Poi 



■ 1 



(4) 



In the case of a plateau, this reduces to Tm/T„.j = 1/Z). These 
results demonstrate that if both the non-subdivided population and 
the isolated deme are in the sequential fixation regime, then 
subdivision significantly accelerates crossing in the best scenario. 
The speedup by subdivision becomes larger (i.e. t„,/t„.5 becomes 
smaller) when the number of demes D is increased at fixed valley 
depth (5 and fixed deme size N (or fixed total population size 
Af = ND). Besides, for Z>» 1, the ratio in Eq. 4 decreases when S is 
increased at fixed N and D: the highest speedups are obtained for 
the deepest valleys. However, as S is increased, the non-subdivided 
population wiU eventually enter the tunneling regime (see Fig. IC). 

Let us now consider the alternative case, where the non- 
subdivided population is in the tunneling regime, while the isolated 
demes are in the sequential fixation regime. In this case, 
Ze = \/(NDpdpm), where pm is the fixation probability of a '1'- 
mutant in an isolated deme (see Erp 1), while T„s = \/(NDpdq), 
where q is the probability that a 'I'-mutant is "successful" in the 
tunneling process, i.e. that its hneage includes a '2'-mutant that 
fixes in the non-subdivided population [28]. Hence, in the best 
scenario, where Tm*Tc, we obtain 



Poi 



(5) 



Since q is independent from population size [28], it also 
represents the probability of successful tunneling in an isolated 
deme. For isolated demes in the sequential fixation regime, q<pm 
by definition [28]. Hence, Eq. 5 entails Xmhm < 1- Thus, speedups 
always exist in the best scenario, provided that the necessary 
condition that isolated demes cross the plateau by sequential 
fixation is satisfied. In the case of a fitness plateau, q = ^/iis [28], 
while poi = 1 /N. Hence, Eq. 5 yields 
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(6) 



nd D 



(10) 



In the other extreme case of a sufTiciently deep valley that 
satisfies &»2^fjis, we have q = fis/S [28]. Using the condition 
(5«1, Eq. 1 yields poi =d/{e^^ — I). Hence, Eq. 5 gives 



= f.is- 



'-1 



(7) 



Interestingly, these expressions oiXmhns are independent of Z) 
at fixed N. This stands into contrast with the regime discussed 
above where the non-subdivided population is in the sequential 
fixation regime. At fixed A'^, the ratio Xmhm expressed in Eq. 7 is 
minimal for 



1.594 



(8) 



The minimum of tmhns, corresponding to the largest speedup 
by subdivision, is obtained for this value of 5: 



: 1.544 A^V- 



(9) 



The small values of mutation probabilities in nature ensure that 
the values of t^/ Tbj in Eqs. 6 and 9 can be very small. 

Conditions for subdivision to maximally accelerate valley 
or plateau crossing 

The previous section was dedicated to the study of the best 
possible scenario, where the valley or plateau crossing time of 
the whole metapopulation is dominated by that, Tc, of the 
champion deme in the metapopulation (i.e. the one that crosses 
fastest). We now determine analytically the conditions under 
which this best possible scenario is attained. For this, we focus on 
migration rates much smaller than division/death rates, 2m<Kd, 
such that fixation or extinction of a mutant lineage in a deme is not 
perturbed by migration. In addition, we assume that isolated 
demes are in the sequential fixation regime, since we showed 
above that it is a necessary- condition for subdivision to significantly 
accelerate crossing, and that it is a sufficient condition for 
subdivision to accelerate crossing in the best scenario. 

In a nutshell, migration must be rare enough for demes to 
evolve quasi-independently, but frequent enough to spread the 
beneficial mutation rapidly. The analytical results below allow for 
predicting the range of migration rates such that subdivision 
maximally accelerates valley or plateau crossing. 

First condition: Quasi-independence. Migration must be 
rare enough for demes to remain shielded from migration while they 
harbor the intermediate mutation. Hence, the average time for a 
deme of 'I'-mutants to fix the beneficial mutation '2', which reads 
Zi2 = l/ri2=l/iNfidpi2), must be smaller than the average 
extinction time, t^, for a deme of T '-mutants to be wiped out by 
migration from other demes with genotype '0'. The total rate of 
migration events in the metapopulation is DNm, so tg =n^,/(DNm), 
where «<. is the average total number of migration events required for 
the '1 '-mutants to go extinct. The first condition, T12 < thus yields 



Let us now estimate He- If one deme has fixed g(-notype '1' while 
all the others still have genotype '0', the probability that a 
migration event involves the mutant deme is Pr = 2/D. Following 
such a "relevant" migration e\'fnt, extinction of the mutant ('1') 
lineage occurs if the '0' migrant fLxes in the '1' deme whUe the '1' 
migrant does not fix in the '0' deme: this occurs with probability 
^10(1 ~Poi)- Conversely, the number of mutant demes increases to 
two with probability ^01 (1 ~Pw)! and otherwise remains constant. 
For Nd»\, using also d«l, we have poi~de^^^«pioxd (see 
Eq. 1). Hence, migration-induced increases in the number of 
mutant demes can be neglected, and we obtain 



1 

PrPlo(l-Poi) 



D 

Yd 



(11) 



In Methods, Sec. 3, we derive the general expression of Tig, 
which does not require Ndy^X, using finite Markov chain theory 
[25]. Note that this general expression is important because 
subdivision generically most accelerates valley crossing for A'^^ « 1 
(see Eq. 8). 

Second condition: Rapid spreading. Migration must be 
frequent enough for the average spreading time of the final 
mutation through the whole metapopulation to be shorter than the 
valley or plateau crossing time T(.~l/(ffoi) by the champion 

deme. Let w.j be the average number of migration events required 
for the final beneficial mutants (with genotype '2') to spread by 
migration, once the champion deme has fixed genotype '2'. Then, 
we can write ts = ns/(DNm), and the second condition reads 



rtsPoi < 



(12) 



Let us now estimate «j, starting from a state where the 
champion deme has fixed genotype '2', whUe all others still contain 
genotype '0'. (Note that some demes may have genotype '1', but 
this is rare since fixation of mutation '1' is the slowest step. 
Moreover, this would not change the spreading time for a plateau 
and would shorten it for a valley.) Let us focus on the regime 
where 5« 1 but A^5» 1, such that mutation '2' is substantially, but 
not overwhelmingly, beneficial [28]. As in the above discussion 
about Tif., we then obtain P20«P02- Thus, it is possible to neglect 
any migration-induced decrease in the number of demes with 
genotype '2', which we denote by /. The probability that a 
migration step exchanges individuals with different genotypes is 
Pi = 2i(D — i)/[D(D—\)], and the probability that such a relevant 
migration step increases ; by one is po2~s. Hence, we obtain 



^' 1 _7)-l-^'l 



DlogD 



(13) 



where the last expr(-ssion is obtained for D»l. In Methods, Sec. 3, 
we use finite Markov chain theory to derive the general analytical 
expression for ris, which does not require Ns»l. 

Combination of the two conditions. Together, Eqs. 10 and 
12 yield the interval of m/ fid over which subdivision maximally 
accelerates valley or plateau crossing. For 
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(14) 



we expect the valley or plateau crossing time T„, of the whole 
metapopulation to be dominated by that of the champion deme: 
'^mhid~^/ D, where TjrfWrQj' is the average crossing time for an 
isolated deme, and in the best scenario, im ic 

In the regime where Ns,N&»\ and .?,^«1, we can use the 
simple expressions of «e and given in Eqs. 11 and 13, which 
yields 



be 



-N5 



^. W 1 

DlogD« — «- 
2 



(15) 



The ratio, R, of the upper to lower bound in Eq. 15 reads 
1 s 



R = 



21) log 



D d\ S) ■ 



(16) 



This ratio increases exponentially with A'^ (this dependence on A'^ 
comes from that of poi)- This entails that, in this regime, the 
interval of m/(p.d) where subdivision most accelerates crossing 
becomes wider as N increases. However, the width of this interval 
is limited by the fact that isolated demes have to be in the 
sequential fixation regime (see Discussion). While the expressions 
of the interval bounds in Er|. 15 are more illuminating and easier 
to derive than the general ones, the latter, given in Methods, Sec. 
3, actually play important roles since the highest speedups of valley 
crossing gained by subdivision are genericaUy obtained for N5 x 1 
(see Eq. 8). 

Case of the fitness plateau. We have obtained an explicit 
expression of the interval of m/fid over which subdivision 

maximally accelerates valley crossing in the case of a relatively 
deep fitness valley where N5 » 1 while (5 « 1 . In the opposite limit 
of a fitness plateau (^ = 0), retaining the assumptions Ns»l and 
i« 1, Eq. 14 can also be simplified. For this, we use the expression 
otrie obtained in Eq. 35 of Methods, Sec. 3, and the expression of 
«.v in Eq. 13, and we note that, since mutation '1' is neutral, 
Poi = \j N and p\2=P02~s. Eq. 14 then becomes: 



(17) 



where we have used iV» 1 and D» 1. The ratio, R, of the upper to 
lower bound in Eq. 1 7 reads 



2„2 



R = 



2D 



(18) 



This simple expression of -R demonstrates that the range of wi/ (pd) 
over which subdivision maximally accelerates plateau crossing 
increases as the deme size N becomes larger, and that this range is 
quite wide as long as the number of demes satisfies D«{Ns)^, 
which is a realistic condition (recall that we are in the regime 
Ns»l). 

Simulation results 

We now present numerical simulations of the evolutionary 
dynamics described above, which enable us to test our analytical 



predictions, and to gain additional insight in the process 
beyond the optimal scenario. Our simulations are based on a 
Gillespie- algorithm [48,49], and described in detail in Methods, 
Sec. 1. 

Let us first focus on the example presented in Fig. ID, which 
shows an example plot of as a function of the ratio of migration 
to mutation rates, m/(jid), obtained through our simulations when 
varying only the migration rate. With the parameter values used in 
this figure, the interx^al of Eq. 14 is 5.8 x 10^^«ff?/(/i<i)«21. 
Note that here, and in the following examples, we use the general 
expressions of ris and given in Methods, Sec. 3, to compute the 
interval of Eq. 14. Fig. ID features a minimum right at the center 
of this theoretically predicted optimal interval. Moreover, this 
minimum corresponds to T,„ =(5.02 + 0.14) x 10', while 
i^W = (3.28 + 0.10) X 10^: hence, the metapopulation crosses the 
valley on average 6.54 times faster than an isolated deme. This is 
very close to the limit of the best possible scenario, where the 
metapopulation would cross 7 times faster than an isolated deme 
(since D = 7 here). This example illustrates that speedups tend 
towards those predicted in the best scenario, when the interval in 
Eq. 1 4 is sulficiently wide (here the ratio between its upper and its 
lower bound is 359). Besides, t„,, = (1.74 ±0.05) x lO*" here: 
comparing it to the above-mentioned value oft,,, yields a 3.47-fold 
speedup of valley crossing by subdivision. The simulation results in 
Fig. ID also show that significant (albeit smaller) speedups exist 
beyond the optimal parameter window. 

Fig. 2 shows heatmaps of the valley crossing time of a 
metapopulation as a function of the migration-to-mutation rate 
ratio, m/(jid) (varied by varying w), and of the fitness valley depth, 
5. Fig. 2A shows that the optimal interval of Eq. 14 (solid lines) 
describes well the region where the ratio Xmhid of the crossing 
time of the metapopulation to that of an isolated deme is smallest 
and tends to the best-scenario limit \/ D. For migration rates lower 
than those in this interval, the ratio im/iid increases when m 
decreases. This can be understood qualitatively by noting that if 
m = Q, is determined by the valley crossing time of the slowest 
among the independent demes. In the opposite case of migration 
rates larger than those in the optimal interval, T„, increases with m, 
and it tends to the non-subdivided case, t„s, at high values of m, as 
expected. Above a threshold value of <5 (dashed line), t„s becomes 
smaller than Ty, in which case large values of m, such that t„ 
tends to t„j, give a low Xmhid (see Fig. 2A). 

Fig. 2B plots the ratio t„,/t„.5 of the crossing time of the 
metapopulation to that of the non-subdivided population, which 
directly yields the speedup obtained by subdividing a population. 
It shows that, for the parameter values chosen, subdivision 
accelerates valley crossing over a large range of vEiUey depths and 
migration rates, extending far beyond the optimal range given by 
Eq. 14, and that the metapopulation can cross valleys orders of 
magnitude faster than a single large population. In addition, above 
a second, larger threshold value of 8 (dotted line in Fig. 2), isolated 
demes enter the tunneling regime [28]: Fig. 2B shows that 
sulficiently above this threshold, the metapopulation no longer 
crosses the valley faster than the non-subdivided population, as 
predicted above. While having isolated demes in the sequential 
fixation regime is a necessary condition to obtain significant 
speedups by sul)di\ision, the non-subdi\idcd j)opulatioii is not 
required to be in the sequential fixation regime (see above, and 
Fig. IC-D). The value of 5 above which the non-subdivided 
population enters the tunneling regime is indicated by a dash- 
dotted line in Fig. 2: significant speedups are obtained both below 
and above this line. The highest speedups are actually obtained 
above it, i.e. when the non-subdivided population is in the 
tunneling regime. With the parameter values used, Eq. 8 predicts a 
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Figure 2. Effect of subdivision on valley crossing time for 
various migration rates and valley depths. A. Heatmap of the 
ratio t„,/t„/ of the average valley crossing time r,„ of a metapopulation 
with D= 10 and A^ = 50 to that r,,/ of an isolated deme with K = 50, as a 
function of valley depth <5 and migration-to-mutation rate ratio m/(^id), 
in logarithmic scale. All numerical results are averaged over 100 
simulation runs, and the heatmap is interpolated. Solid lines: bounds of 
the interval in Eq. 14. Dashed line: value of S above which a non- 
subdivided population crosses the valley faster than an isolated deme. 
Dotted line: value of d above which an isolated deme is in the tunneling 
regime. Dash-dotted line: value of d above which the non-subdivided 
population is in the tunneling regime. Parameter values: a' = 0.1, 
^1 = 5 X 10^^, .5 = 0.3; (5 and m are varied. B. Similar heatmap for the ratio 
Tm/x,!., of the average valley crossing time t„, of a metapopulation to 
that T„., of a non-subdivided population (with if = 500). Solid line: 
predicted value of d, from Eq. 8, for which the largest speedup by 
subdivision is expected. 
doi:1 0.1 371/journal.pcbi.1 003778.g002 



minimum of t,„/t,„ for 0.035 (solid line in Fig. 2B), which 
agrees very well with the results of our numerical simulations. 
(Note that this value of 8 satisfies Sy>2y^, and is such that the 
non-subdivided population is in the tunneling regime. These 
conditions were used in our derivation of Eq. 8.) 

Discussion 

Limits on the parameter range where subdivision 
maximally accelerates crossing 

In the Results section, we have shown that having isolated 
demes in the sequential frxation regime is a necessary condition for 
subdivision to significantly accelerate crossing. This requirement 
limits the interval of the ratio m/(j.id) over which the highest 
speedups by subdivision are obtained. The extent of this interval 
can be characterized by the ratio, R, of the upper to lower bound 
in Eq. 14. Let us express the bound on R imposed by the 
requirement of sequential fixation in isolated demes. 

If 2y/Is«(5«l, the threshold value A^x below which an 
isolated deme is in the sequential fixation regime satisfies 
e^^^xS^/dis) [28]. Let us also assume that NS»l, and that 
s«l while iVj'»l, to be in the domain of validity of Eqs. 15 and 
16. Combining the condition N <Nx with the expression of in 
Eq. 16 yields 



5 (\ ^ 
2fiD\ogD V ^6 



(19) 



For plateaus, isolated demes are in the sequential frxation 
regime if their size N is smaller than Nx = l/y^ [28]. In the 
regime of validity of Eqs. 17 and 18 (.?« 1 while Ns»l, and 1, 
D»l), this condition can be combined with Eq. 18, which yields 



R< 



2fiD ' 



(20) 



Both Eq. 19 and Eq. 20 show that increasing the number D of 
demes decreases the range where the highest speedup by 
subdivision is reached. This is because having more subpopula- 
tions makes the spreading of the beneficial mutation slower. In 
addition, we find that the bound on R is proportional to \/fi. 
Hence, despite this bound, the interval where subdivision most 
accelerates plateau crossing can span several orders of magnitude, 
given the small values of the actual mutation probabilities ^ in 
nature. 

Effect of varying the degree of subdivision of a 
metapopulation 

An interesting question raised by our results regards the optimal 
degree of subdivision. Given a certain total metapopulation size, 
into how many demes should it be subdivided in order to obtain 
the highest speedup possible? We first attack this question using 
our analytical results, and then we present simulation results, 
which allow for going beyond the best scenario and its associated 
parameter window. 

Let us consider a metapopulation of given total size J\f = ND. 
Our analytical results show that increasing subdivision, i.e. 
increasing the number D of subpopulations at constant A/", leads 
to stronger speedups of valley crossing (see Eqs. 4 and 7, with 
N=J\f/D). However, Eqs. 16 and 18, and the previous 
paragraph, show that when D is increased, the parameter range 
where the speedup by subdivision tends to the best-scenario value 
becomes smaller and smaller. Eventually, this parameter range 
ceases to exist altogether: this occurs when R becomes of order 1 
and below. This sheds light on an interesting trade-off in the 
degree of subdivision D, between the magnitude of the optimal 
speedup gained by subdivision and the width of the parameter 
range over which the actual speedup is close to this optimal value. 
This effect can be observed qualitatively in Fig. 3A, where the 
valley crossing time T,„ of a metapopulation with fixed total size is 
shown versus the migration-to-mutation rate ratio, m/{^id), for 
different values off: when D is increased, the minimum becomes 
deeper but less broad. 

In addition, Eqs. 15 and 17 show that when D is increased, the 
lower bound of the interval where the speedup by subdivision 
tends to the best-scenario value decreases, as D log D for plateaus 
(Eq. 17) and even more rapidly for deep valleys (Eq. 15). 
Qualitatively, this is because spreading of the beneficial mutation 
gets longer when D increases. Conversely, the upper bound of this 
parameter range is independent of D for deep valleys (Eq. 1 5), and 
grows only logarithmically with D for plateaus (Eq. 17). Hence, 
when D is increased, the center of the interval where the actual 
speedup is close to the optimal value shifts towards higher 
migration rates. This effect, which can be observed in Fig. 3A, is 
studied more precisely in Fig. 3B: at fixed migration rate m, the 
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Figure 3. Varying the degree of subdivision of a metapopulation. A. Valley crossing time t„, of a metapopulation with total carrying capacity 
DK = 2500, versus migration-to-mutation rate ratio m/{pid), for four different numbers D of demes. Dots are simulation results, averaged over 1000 
runs for each value of m/{iid) (500 runs for a few points far from the minima); error bars represent 95% CI. Vertical lines represent the limits of the 
interval of m/{fid) in Eq. 14 in each case, except for £)= 125, where this interval does not exist. Black horizontal line: plateau crossing time for a non- 
subdivided population with .S^ = 2500 for the same parameter values, averaged over 1000 runs; shaded regions: 95% CI. Dashed line: corresponding 
theoretical prediction from Ref. [28]. Parameter values: d = 0.l, ^i = 8x 10^', .? = 0.3 and S = 6x 10^^ (same as in Fig. 1C-D); m is varied. B. Valley 
crossing time t,„ of a metapopulation with total carrying capacity DK = 2500, versus the number D of demes, for m = 10^^ (i.e. m/{i.id) = 12.5). Dots 
are simulation results, averaged over 1000 runs for each value of D; error bars represent 95% CI. Parameter values: same as in A. C. Valley crossing 
time T,nin, minimized over m for each value of D, of a metapopulation with total carrying capacity DK = 2500, versus the number D of demes. For 
each value of D, the valley crossing time of the metapopulation was computed for several values of m, different by factors of 10° or 10° in the 
vicinity of the minimum (see A): T„,in corresponds to the smallest value obtained in this process. Results obtained for the actual metapopulation (blue) 
are compared to the best-scenario limit (red) where r,nin = t„//A calculated using the value of t„/ obtained from our simulations. Dots are simulation 
results, averaged over 1000 runs for each value of D; error bars represent 95% CI. Dashed line: value of D such that .R= 100. Dotted line: value of D 
above which the deleterious mutation is effectively neutral in the isolated demes. Solid line: value of D such that R=l. Parameter values: same as in 
A and B. 

doi:1 0.1 371 /journal.pcbi.1 003778.g003 



crossing time T„, of a metapopulation exhibits a minimum at an 
intermediate value of D. Indeed, the crossing time of the 
metapopulation first decreases when D is increased because the 
minimum crossing time then decreases. But beyond a certain value 
of D, the migration rate that yields the highest speedup becomes 
larger than the fixed migration rate m, so T,„ increases when D is 
increased further. 

Next, we study the dependence on D of the valley crossing time 
Tmin minimized over in for each D, again for a metapopulation 
with fixed total size A/" = ND. For values of D small enough for the 
interval in Eq. 14 to be broad, we expect Tn,in to be close to the 
optimal scenario value T,rf/-D. But, as discussed above, as D 
increases, this interval wiU become smaller and then vanish. In 
such a regime, our analytical results are no longer sufficient to 
predict the dependence of T„iin on D, but our simulations can 
provide additional insight. Fig. 3C shows that, while _R> 100 (left 
of the dashed line), Tnii„ is close to the best-scenario value. When D 
is increased beyond this point, Tmin decreases slower than the best- 
scenario value. Indeed, the interval in Eq, 14 is no longer wide 
enough for the best-scenario limit to be approached. Note also that 
when demes become small enough, verifying N=J\f/D«l/5 
(right of the dotted line in Fig. 3C), mutation '1' becomes 
effectively neutral in individual demes, as ^oi tends to l/N (see 
Eq. 1). For even higher values of D, Tmin is observed to saturate 
rather than exhibiting a unique minimum. Interestingly, this 
occurs for D such that the interval in Eq. 14 fuUy vanishes (i.e. 
when R passes below 1, right of the solid line on Fig. 3C). While 
we do not have rigorous proof of the generic existence of this 
saturation, we have explored this point for other parameters, and 
found similar behavior (data not shown). Importantiy, this 
indicates that there is a whole class of nearly optimal population 
structures. 



Extension to weakly beneficial intermediates 

Our work has focused on fitness valleys (<5>0), such diat 
mutation ' 1 ' is deleterious, and on fitness plateaus (5 = 0), such that 
mutation '1' is neutral. For \S\ < max (y^.l/A^), mutation '1' is 
effectively neutral, as far as valley crossing is concerned, in a 
population with individuals [28]. (This condition holds both in 
the sequential fixation regime and in the tunneling regime.) This 
implies that our arguments and our results obtained in the case of 
the fitness plateau also hold for weakly beneficial intermediates. 
This point is illustrated in Fig. 4A. 

Extension to a population coupled to small island 
populations 

Thus far, we focused on demes of equal size for simplicity, but 
demes of different sizes are relevant in practice. As a step toward 
more general populations structures, we now consider a popula- 
tion connected by migration to S smaller satellite populations of 
identical size, assumed to be in the sequential fixation regime. We 
only allow migration between the large population and each of the 
smaller islands, and the total migration rate is denoted by M. The 
small island affected by migration is chosen randomly at each 
migration event. It is straightforward to adapt our work to this case 
(see Methods, Sec. 4). We obtain an interval of M/{jid) over 
which the crossing time for the large population is dominated by 
the crossing time of the champion island. This is corroborated by 
our simulations (see Fig. 4B). 

Realistic parameter values 

Let us consider the example ot Escherichia coli, for which the 
mutation probability per base pair per division is/i«8.9xlO^" 
[50]. In order to gain a speedup of crossing by subdivision, we 
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Figure 4. Extension to effectively neutral intermediates and to 
a large population connected to smaller islands. A. Valley 
crossing time t,„ of a metapopulation composed of /)=10 demes with 
K = 130, versus migration-to-mutation rate ratio m/(iid), in logarithmic 
scale, for three values of <5 in the effectively neutral regime. Dots are 
simulation results, averaged over 100 runs for each value of m; error 
bars represent 95% CI. Black vertical lines represent the limits of the 
interval of m/(nd) in Eq. 14. Black horizontal line: plateau crossing time 
for an isolated deme with if =130 for the same parameter values, 
averaged over 100 runs; shaded regions: 95% CI. Dashed line: 
corresponding theoretical prediction from Ref. [28]. Note that the 
plateau crossing time of the non-subdivided population with ^ = 1300 
is indistinguishable from that of the isolated deme here (as both are in 
the sequential fixation regime). Parameter values: d = OA, /( = 5 x 10-', 
5 = 0.5; m is varied. B. Valley crossing time i/ of a large population with 
A^ = 500 connected to S smaller islands with K' = 50, versus migration- 
to-mutation rate ratio M/did), in logarithmic scale. Dots represent 
simulation results averaged over 100 runs for each value of M, and 
error bars are 95% CI. Vertical lines represent the limits of the interval of 
Ml(pd) in Eq. 49. Blue (resp. red) line: valley crossing time for an 
isolated population with = 50 (resp. K = 500) for the same parameter 
values, averaged over 100 runs; shaded regions: 95% CI. Dashed blue 
(resp. red) lines: corresponding theoretical predictions from Ref [28]. 
For 5=1, the observed minimum t/ satisfies t/5;t„, where t„ is the 
average valley crossing time of an isolated island. For 5=10, the 
observed minimum satisfies t/s;t„/10. Parameter values: a' = 0.1, 
/i = 4x 10-^ 5 = 0.25 and ^ = 0.1; M is varied. 
doi:1 0.1 371/journal.pcbi.1 003778.g004 



require demes to be in the sequential fixation regime. For plateaus, 
this condition reads N <\/ Let us consider deme sizes such 
that this condition is satisfied. 

First, let us choose A'^=5x IC*, which is within the smallest 
range of sizes used in current evolution experiments. For instance, 
it is the number of bacteria transferred at each dilution step for 
small populations in [27]. For this value of A^, all plateaus with 
5<4.5 are in the sequential fixation regime (from the condition 



N < \/ ,JJIs). Let us also consider D= 100, since 96-weU plates are 
often used in these experiments [27,51]. This yields a total 
population size of 5 x 10^ individuals, which is in the tunneling 
regime for all plateaus with 5>4.5 x IQ—*. For s= 10-^, isolated 
demes are in the sequential fixation regime for 0< (5< 2.2 x lO-''. 
(Subdivision cannot significantiy accelerate crossing for deeper 
valleys since isolated demes are then in tunnehng, but those valleys 
take longer to cross than shallow ones and are thus probably less 
often crossed in practice.) The ratio R of the bounds of the interval 
in Eq. 14 satisfies R>'i25 throughout this range of valleys, with 

« 10' for the plateau and lO'' for the deepest valleys in the 
range. Thus, actual speedups will approach the best-scenario one, 
and significant speedups will exist in a wide parameter window. 
Eq. 8 predicts that the highest speedup is obtained for 
(5«3.2xlO-', and Eq. 9 then yields a speedup factor by 
subdivision of t„.v/t,„ « T„s/Tf. «2.9 x 10^. (Using instead the full 
expression of obtained from Eq. 23 (see Methods, Sec. 2) yields 
2.7 x 10^, i.e. a correction of 7%.) Moreover, for all valleys with 
(5<3.2xlO-', the best-scenario speedup ranges from 18 to 
2.7 x 10^. Thus, subdivision significantly accelerates crossing for 
this entire class of valleys. 

It should be noted that the timescales obtained in this example 
are long compared to experimental ones. For instance, for the 
plateau, t^. corresponds to 1.3 x 10* divisions while T,„ is 2.4 x 10' 
divisions. However, Tf. can become smaller if the number of 
subpopulations D is increased, as discussed in our previous section. 
Besides, we have chosen to focus on standard Escherichia coli for 
simplicity. Organisms with a higher mutation rate, e.g. viruses 
such as HIV, or mutator strains, would have much shorter 
timescales, but smaller subpopulations would then be required for 
demes to be in the sequential fixation regime. 

Our example thus far focused on a small but realistic deme size, 
A^ = 5 X 10'*. Experimentally more frequent values of are in the 
range 5 x 10^ - 10' [27,51]. Increasing N at fixed n decreases the 
range of s for which demes are in the sequential fixation regime. 
For a plateau, this condition reads s< l/{fiN^). For N = 5 x 10^, 
this yields 5<4.5xl0-^, and for A^ = 5xlO*, this yields 
.y<4.5 X lO-"*. Hence, the range of plateaus (and similarly, of 
valleys) for which subdivision accelerates crossing becomes more 
restricted when N is increased. Nevertheless, if these increasingly 
stringent conditions on s are satisfied, significant speedups by 
subdivision are still expected. Indeed, Eq. 9 shows that the smallest 
value of the ratio Tc/t„j is proportional to N^s, so if one increases 
while decreasing s as 1/iV^, the maximal speedup by 
subdivision will remain unchanged. 

In this work, we have considered the crossing of one particular 
valley or plateau corresponding to a specific pair of two mutations. 
Given the complexity and high dimensionality of actual fitness 
landscapes, there may be a large number of parallel valleys or 
plateaus, so that one of these could be crossed quite frequently 
even though the crossing time for a single valley or plateau 
remains large. Our work shows that, under specific conditions, 
subdivision can significantly accelerate crossing for whole classes of 
valleys and plateaus. Furthermore, in a generic, high-dimensional 
fitness landscape that contains both valleys and/or plateaus and 
uphill paths, subdivision can provide an additional effect: it 
'^shields" some demes in the metapopulation from adaptation via 
the uphill paths, leaving them time to explore valley-crossing paths 
that may be better in the longer term. While this effect is outside 
the scope of the present paper, it could lead to additional 
advantages of subdivision in evolution on rugged fitness land- 
scapes. 
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Conclusion 

Our study of a generic and minimal model of population 
subdivision with migration demonstrates that subdividing a 
population into demes connected by migration can significantiy 
accelerate the crossing of fitness plateaus and valleys, without the 
need for additional ingredients. We have derived quantitative 
conditions on the various parameters for subdivision to accelerate 
crossing, and for the resulting speedup to be maximal. In 
particular, isolated demes have to be in the sequential frxation 
regime for a significant speedup to occur. This condition is quite 
strong, but provided that it is met, significant speedups can be 
obtained in a wide range of migration rates, with the fastest deme 
driving the crossing of the whole metapopulation in the best 
scenario. We have derived the interval of migration rates for which 
this best scenario is reached. In addition, we have shown that 
increasing the degree of subdivision of a population enables higher 
speedups to be reached, but that this effect can saturate. 

Our quantitative assessment of the conditions under which 
subdivision significantly speeds up valley or plateau crossing can 
aid in optimally designing future experiments, enabling one to 
choose the sizes and the number of demes, as well as the migration 
rates, such that subdivision can accelerate valley and plateau 
crossing. 

Further directions include investigating the evolution of a 
metapopulation with a distribution of deme sizes on a more 
general rugged landscape, as well as assessing the impact of specific 
geographic structure. Our work could also be extended to sexual 
populations, where recombination plays an important role in 
valley or plateau crossing [52]. The interplay between recombi- 
nation and subdivision, which respectively alleviate and exacerbate 
clonal interference, would be interesting to study. 

Methods 

1 Simulation methods 

Our simulations are based on a Gillespie algorithm [48,49] that 
we coded in the C language. Here we will describe our algorithm 
for the case of a metapopulation of D demes of identical size, 
which is the primary situation discussed in our work. In our 
simulations, each deme has a fixed carrying capacity K-we discuss 
this choice further in this section. 

1.1 Algorithm. A number of different events occur in our 
simulations, each with an independent rate: 

• Each individual divides at rate fg(l—Ni/K), where fg is the 
fitness associated with the genotype g e {0,1,2} of the 
individual, and A^, is the current total number of individuals 
in the deme / e [i,D\ to which the individual belongs. This 
corresponds to logistic growth. 

• If a dividing cell has g<2, upon division, its offspring (i.e., one 
of the two individuals resulting from the division) mutates with 
probability /i, to have genotype g+l instead of g. 

• Each individual dies at rate d. Hence, at steady-state, 
NiXK{l—d/fi), where / is the average fitness of deme /. In 
practice, we choose d = OA, and fitnesses of order one, thus 
NiK0.9K. 

• Migration occurs at total rate m i -^i- Two different demes 
are chosen at random, an individual is chosen at random from 
each of these two demes, and the two individuals are 
exchanged. There is no geographic structure in our model, 
i.e. exchange between any two demes is equally likely. 

In practice, the number of individuals with each genotype in 
each deme is stored, as well as the corresponding division rate. 



This data fuUy describes the state of the metapopulation, and 
allows determination of the rates of all events. For each event in 
the simulation, the following steps are performed: 

• A timestep dt is drawn from an exponential distribution with 
rate equal to the total rate R, of events (i.e., the sum of all 
rates), and time is increased from its previous value, /, to t + dt. 
In other words, the next event occurs at time t + dt. 

• The event that occurs at t + dt is chosen randomly, in such a 
way that the probabihty of an event with rate r is equal to 
r/R,: either a cell divides, or a cell dies, or a migration event 
occurs. 

• The event is performed, and the relevant data is updated. 
Since we store the number of individuals with each genotype in 
each deme, only one or two of these numbers need to be 
updated at each step. In addition, the division rates of the 
affected deme i must be updated upon division and death 
because A,- is modified. Note, however, that this represents 
only three numbers at most (one for each genotype). 

The advantage of the Gillespie algorithm is that it is exact, and 
does not involve any artificial discretization of time. 

1.2 Working at fixed carrying capacity. In our simulations, 
demes have a fixed carrying capacity, and the number of 
individuals per deme fluctuates weakly around its equilibrium 
value. This approach, also used in e.g. [23], has the advantage of 
realism. Alternatively, we could impose a constant number of 
individuals per deme. 

(i) First, we could choose a dividing individual in the whole 
metapopulation with probability proportional to its fitness, 
and simultaneously suppress another individual, chosen at 
random in the same deme. However, in this case, individuals 
in demes of higher fitness would exhibit shorter lifespans, 
which is not realistic and may introduce a bias. 

(ii) A second possibility would be to choose a dividing individual 
(according to fitness) in each of the demes, and to 
simultaneously suppress another individual, chosen at 
random, in each deme. However, in this case, unless 
migration events are far less frequent than these collective 
division-death events (i.e., these D division-death events), the 
time interval between them becomes artificially discretized. 
This introduces biases unless the total migration rate mDN is 
much smaller than Nd, i.e. unless m<iid/D. 

Consequently, while imposing a constant number of individuals 
is a good simulation approach for a non-subdivided population 
(see e.g. [28]), it tends to introduce biases in the study of 
metapopulations. While we chose to perform simulations with 
fixed carrying capacities in order to avoid any of these biases, we 
checked that, for small enough migration rates, our results are 
completely consistent with simulation scheme (ii) described above. 
This consistency check also demonstrates that it is legitimate to 
compare our simulation results obtained with fixed carrying 
capacities to our analytical work carried out with constant 
population size per deme. 

2 Crossing time of the champion deme 

In this section, we give more details on the calculation of the 
average valley or plateau crossing time Tj. by the champion deme 
amongst D independent ones. We show in the Results section that, 
in the best scenario, the crossing time of the whole metapopulation 
is determined by this time. 
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Zc is the average shortest crossing time of Z) independent demes. 
This minimum crossing time, which we denote by tc, is also called 
the smallest (or first) order statistic of the deme crossing time 
amongst a sample of size D [53]. 

Let us denote by p(i) the probabihty density function of valley 
or plateau crossing time for a single deme, and let us introduce 



V{t)-- 



p(t')dt' (it satisfies 'P{t) = \- C{t) where C(t) is the 



cumulative distribution function of valley or plateau crossing by a 
single deme). The probability that tc is larger than / is equal to the 
probability that the crossing times of each of the D independent 
demes are all larger than t : P{tc>t) = \P(t)Y' . By differentiating 
this expression, one obtains the probability density function Pc(tc) 
of the crossing time tc by the champion deme (see e.g. [53]): 



Pc(.tc) = D[V(tcT-'p{tc). 



(21) 



We now express p{t) exphcidy. Since demes are assumed to be 
in the sequential frxation regime, valley or plateau crossing 
involves two successive steps. The first step, fixation of a 'V- 
mutant, occurs with rate roi, and the second step, frxation of a '2'- 
mutant, occurs with rate r\2 (see the Results section for expressions 
of these rates). The total crossing time is thus a sum of two 
independent exponential random variables, with probability 
density function given by a two-parameter hypoexponential 
distribution [53]: 



pity- 



''""•'^ (e-'oi'- 

''12 -''01 



(22) 



Combining Eqs. 21 and 22, we obtain 

^ri2e"'oi''--foie"''i2'''^ 



Pc(tc) = D 



''12— ''01 



P(.tc), (23) 



with p(tc) given by Eq. 22. tc can then be determined for any 
value of the parameters by computing the average value of tc over 
this distribution. 

Since mutation '1' is deleterious or neutral while mutation '2' is 
beneficial, the first step of valley crossing is much longer than the 
second one over a broad range of parameter values. In this case, 
we can approximate p(t) with a simple exponential distribution, 



p{t) = me 



Eq. 2 1 then yields 



Pc(tc) = Drme 



(24) 



(25) 



i.e. tc is distributed exponentially with rate Dyqi . In this case, we 
simply have Tf.« l/(Z)roi), which can be written as Zc~Xici/D, 
where Xjd « 1 /roi is the average crossing time for an isolated deme. 
Hence, in this case, on which our analytical discussion focuses, the 
champion deme crosses the valley D times faster on average than 
an isolated deme. 

For this approximation to be valid, the second step of valley 
crossing must be negligible even for the champion deme, i.e., 
Dpoi «pi2- For very large D, the value of will not be as small as 



l/(Z>foi), since the second step will no longer be negligible (see 
[52] for a discussion of similar issues). The crossover to this regime 
can be determined by computing the average of the distribution in 
Eq. 23 and comparing it to l/(Z)foi). 

3 Number of migration events for extinction or spreading 
in a metapopulation 

In our Results section, we have derived an interval of the ratio 
of migration rate to mutation rate over which subdivision most 
reduces valley or plateau crossing time (see Eq. 14). The upper 
bound involves Mp, the average number of migration events 
required for the 'I'-mutants to be wiped out by migration, starting 
from a state where one deme has fixed genotype '1', while all other 
demes have genotype '0'. Similarly, the lower bound involves n^, 
the average number of migration events required for the '2'- 
mutants to spread by migration to the whole metapopulation, 
starting from a state where one deme has fixed genotype ^2', while 
all other demes have genotype '0'. In our Results section, we have 
provided intuitive derivations of the simple expressions of «j and 
He, valid for Ns»l and s«l, NS»l and 5«\ (see Eq. 15). 
However, it is important to derive more general expressions, 
especially since subdivision genericaUy most accelerates vaUey 
crossing in the intermediate regime where N3 x 1 (see Results, Eq. 
8). 

Here, we derive general analytical expressions for and «j, 
both for fitness plateaus and for fitness valleys. These more general 
expressions are those used for numerical calculations of the bounds 
in our examples. Throughout this section, we consider a 
metapopulation of D demes composed of N individuals each, 
and we assume that individual demes are in the sequential frxation 
regime (see Results). 

3.1 A finite Markov chain. In order to determine and M^, 
we study the evolution of the number / e [0,D] of demes that have 
fixed the mutant genotype ('1' for the calculation of «(,; '2' for that 
of Mj), while other demes have genotype '0'. Given that the value of 
i just before a migration step fuUy determines the probabilities of 
the outcomes of this migration step, and given that i = 0 and i = D 
are absorbing states, the number ; evolves according to a finite 
Markov chain, each step being a migration event. We next express 
the transition matrix of this Markov chain. 

The only migration events that can affect i are those that 
exchange individuals from two demes with different genotypes. Let 
us call these migration events "relevant". The probability p'- of a 
migration event being relevant corresponds to the probability that 
this migration affects one of the i mutant populations and one of 
the D — i '0' populations: p- =2i{D — i)/[D{D— 1)]. We only focus 
on the final outcome of a migration event, after frxation or 
extinction of each of the two migrants' lineages has occurred. Let p 
denote the probability that the mutant migrant fixes in the '0' 
deme, and p' the probability that the '0' migrant fixes in the 
mutant deme. As a result of one such relevant migration event: 

• i increases by one with probability p{l—p'), if the migrant 
mutant fixes in the ^0' deme while '0' migrant does not fix in 
the mutant deme. 

• i decreases by one with probabihty />'( 1 —p), in the opposite 
case. 

• Otherwise, / does not change. This happens either if both 
migrants fix (with probability pp') or if no migrant frxes (with 
probability (1 —p)( l -/>'))■ 

These probabilities, multiplied by the probability p'- that a 
migration event is relevant, yield the transition matrix of our finite 
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Markov chain, which is tri-diagonal (or continuant) since each 
migration step can either leave ; constant, or increase or decrease 
it by one: 



2i{D-i) 
D(D-l) 



Pil-P'), 



(26) 



vo = 



( 


Ea-.; Pkj 


1 

1 


{HkJ Pk) 


^{T.k=! Pk] 


\P]Pj^J+l 



(31) 



Pi^i- 1 = i^-P)' (27) y= 1 (^E*:=o j PjPj^j+ 1 



where we have introduced 



(28) 



for (e[l,£) — 1], and Pq^q=Pb^d = 1- We have denoted by Pj^k 
the probability that / varies from / to k as the final outcome of one 
migration event. 

Here, we do not account for independent mutations arising and 
fixing in other demes during the process of spreading (or 
extinction) of the mutant's lineage in the metapopulation. Indeed, 
our aim is to compare the timescales of migration and mutation 
processes, so we treat them separately. Note that, in practice, this 
hypothesis is reasonable if mutations that fix are sufficientiy rarer 
than migration events. We also consider that the time between two 
successive migration events is large enough for fixation to occur in 
the demes affected by migration before the next migration event 
occurs, which is true in the low-migration rate regime that we 
study in our work {2m «d, where 2m is the migration rate per 
individual, while d is the death and division rate per individual). 

ng and can be directiy expressed as the average number of 
steps of the Markov chain necessary to go from the initial state 
i = l to absorption in a particular absorbing state, either i = D or 
/ = 0. Let us present general expressions of these average numbers 
of steps, before using them to obtain explicit expressions of Wj and 

tte. 

3.2 Some results regarding finite Markov chains with tri- 
diagonal probability matrices. We are interested in the 
average number of steps Va until the system reaches each of the 
absorbing states a e {0,D}, starting from the state /= 1: 



D-l 



Va= X^-S/,a. 



(29) 



Pk-- 



k p. . . 



(33) 



3.3 Explicit expression of n^. tie, in fact, corresponds to Vq, 

where p is the probability that a ' 1 '-mutant fixes in a deme of '0' 
individuals (i.e. p=pQi) and p' is the probability that a '0'- 
individual fixes in a deme of 'I'-mutants (i.e. p' =pio)- Hence, it 
can be expressed explicitiy from Eqs. 31, 26, and 27. Since the 
expressions of p and p' depend whether mutation '1' is neutral or 
deleterious, we obtain different expressions for the fitness plateau 
and for the fitness valley. 

Fitness plateau. For a fitness plateau (i.e. a neutral 
intermediate '1'), p=p' = l/N, where N is the number of 
individuals per deme. Hence, 



_ _ 2KD-i)iN-l) 
Pt^t+ 1 - Pi^i- 1 - ' (34) 

which implies that ft = 1 for all A: (see Eq. 33). Thus, Eq. 31 yields 



N-D 



■E 



1 N 



2(N^\)jrij 



DlogD, 



(35) 



where the last expression holds for A'^» 1 and Z)» 1. 
Fitness valley. Eqs. 33 and 26, 27 yield Pi^ = r^, with 



a-p)p' 
ii-p')p' 



(36) 



and Eq. 31 gives: 



where is the average number of steps that the system spends in 



the state i = j before absorption, given that it starts in the state / = 
and finally absorbs in state i = a. It can be expressed as [25] 



I 



2(r- 



-r^)(l-0(l- 



PlPfrlrJjXD-j)- 



(37) 



(30) 



where Sj is the average number of steps the system spends in state 
i=j before absorption in either of the two absorbing states, given 
that it started in state i = l, and nj^a is the probability that the 
system finally absorbs in state i = a if it starts in state / = j. 

Using the explicit expressions given in [25] for Sj and 7iy_a in the 
case of a tri-diagonal probability matrix, we obtain: 



In these expressions, p=pm is the probability of fixation of a 
deleterious '1 '-mutant, with fitness 1 —(5, in a deme where all other 
individuals have genotype '0' and fitness 1. It can be obtained from 
Eq. 1 , as well as the probability p' = pio of the opposite process. 

3.4 Explicit expression of «j. corresponds to Vr>, where 
P=P02 is the probability that a '2'-mutant (with fitness 1-1-5) fixes 
in a deme of '0' individuals (with fitness 1), and p'=P20 is the 
probability that a 'O'-individual fixes in a deme of '2'-mutants. 
Hence, it can be expressed explicitiy from Eqs. 32, 26, and 27, 
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using Eq. 1 to express the fixation probabilities. For rig, there is no 
difierence between the valley and the plateau, since genotype 'I' is 
not involved. 

As above, Eqs. 26, 27 and 33 yield p/^ = i^, with r defined in Eq. 
36. Thus, Eq. 32 gives 



D(D-l) 



V 

_ n'\n ^—z 



(l-rJ){l-r^-J) 



(38) 



3.5 Simplified expressions for deep valleys and for 
plateaus. In our Results section, we have shown that the 
benefit of subdivision is highest when m/(jid) is situated between a 
lower bound, 



L = nsPoi, 



and an upper bound. 



U-- 



D 



(39) 



(40) 



(see Eq. 14), where p\2 denotes the probability of fixation of a 
single mutant with genotype '2' in a background of 'I'-mutants. 
Here we present simplified expressions for Mj and rig, and hence of 
L and U , in particular parameter regimes. 

Throughout this section, we focus on the regime where Ns » 1 
but J« 1, such that mutation '2' is substantially, but not 
overwhelmingly, beneficial [28]. We then have /)20~'56~*«1 
and po2~s (see Eq. 1). To leading (i.e. zeroth) order in /720, we 
obtain from Eq. 38 that 



Z»-l^'l D\ogD 



j=lJ 



(41) 



where the last expression holds for Z>»1. This expression of is 
identical to Eq. 13, which was demonstrated more intuitively in 
the Results section by directly assuming Ns » 1 and s « 1 . 

We now consider the case of a plateau (d = 0) and the case of a 
valley such that N5 » 1 but S«l. We demonstrate that the latter 
case is consistent with the simplified derivations in our Results 
section. 

Fitness plateau. For a fitness plateau, combining Eqs. 35 and 

40, the upper bound U reads 



Ns, „ 
t/w— logZ), 



(42) 



where we have used pn=PQ2 since mutation '1' is neutral, and 
assumed N» 1 and Z)» 1. 

Additionally, Eqs. 41 and 39 can be combined to write the 
lower bound L as 



D 

'Ns 



logD. 



(43) 



again to lowest order in p2o- Here, we have used poi = l/N, since 
in the case of the plateau, mutation ' 1' is neutral. This expression 
too holds for A'^» 1 and D»l. 

Combining Eqs. 42 and 43 yields Eq. 17. 

Fitness valley. Next we focus on valleys such that d«l but 
Nd»l. (Note that, in the opposite limit Nd«l, mutation '1' is 



effectively neutral, and the above discussion regarding the fitness 
plateau applies.) Then, poi xSe^^^ «l (see Eq. 1). To lowest (i.e. 
zeroth) order in poi, Eq. 37 becomes 



D 

2pio 



D 

'Yd' 



(44) 



where we have used the approximation P\qK&, which holds for 
8« \ and Af(5»l. This expression of coincides with Eq. 11, 
which is obtained in the Results section through a more intuitive 
argument that directly assumes 5« 1 and Nd'»\. Hence, from Eq. 
40, the upper bound U is 



U', 



Pu 
' lb 



(45) 



where we used the conditions 5«1, Nby>\, J«l and Nsy>\ to 
simplify the expression oi p\2. 

Meanwhile, from Eq. 39 and 41, the lower bound L takes the 
form 



L = ^i)log/). 

P^i 



be 



-m 



DIogD, 



(46) 



where, again, we used the conditions <5«1, A'^<5»1, s«l and 
Ns»l to simplify the expressions of/>oi and po2- 
Combining Eqs. 46 and 45 yields Eq. 15. 

4 A population connected by migration to smaller 
population islands 

Let us consider a population of Af individuals connected by 
migration to S smaller population islands with N <J\f individuals 
each. These islands of identical size are assumed to be in the 
sequential frxation regime. For the sake of simplicity, we consider 
that migration only occurs between the large population and the 
islands: a migration step is a random exchange of two individuals 
between the large population and one of the islands (chosen at 
random at each migration event), and the total migration rate is 
denoted by M. Here, we focus on the valley or plateau crossing 
time of the large population. We demonstrate that the evolution of 
a large population can be driven by that of satellite islands. 

In the optimal case, the crossing time of the large population is 
determined by that of the champion island, i.e., that which crosses 
the fitness valley or plateau fastest. We now determine the 
conditions under which this optimum is achieved, focusing on 
migration rates much smaller than division/death rates, 
M« min (dNS,dJ\f), such that fixation or extinction of a mutant 
lineage in either the large population or an island is not 
significantly perturbed by migration. Again, migration should be 
rare enough for islands to remain effectively shielded from 
migration events while they have fixed the intermediate mutation, 
until the final beneficial mutation arises. Second, migration should 
also be frequent enough for the spreading time of the final 
beneficial mutation from the champion island to the large 
population to be negligible with respect to the crossing time of 
the champion island. These two criteria again provide upper and 
lower bounds on M/{p.d). 

The average time Ti2 = I /(Nfidp 12) (with pu from Eq. 1) 
required for an island of ' 1 '-mutants to fix the beneficial mutation 
'2' must be smaller than the average time, tg, for an island of '1'- 
mutants to be wiped out by migration from the large population, 
which stiU exhibits genotype '0'. The rate of migration events 
between the island of '1 '-mutants and the large population is 



PLOS Computational Biology | www.ploscompbiol.org 



13 



August 2014 I Volume 10 | Issue 8 | el 003778 



Population Subdivision and Rugged Landscapes 



M/S. Hence, te = S/(Mpio), where pio is the probability of 
fixation of the lineage of a single migrant with genotype '0' in an 
island where all other individuals are 'I'-mutants: for valleys, it is 
given by Eq. 1, while for plateaus, it is equal to 1/A'^. The first 
condition, T12 < ?e, thus yields 



M NSpn 



(47) 



The second condition is that the average spreading time, 1^, for 
the final beneficial mutation to fix in a large population after it has 
fixed in the champion island, must be smaller than the average 
valley or plateau crossing time, T^, of the champion island. Similar 
to te previously, we obtain ts = S/(MpQ2), where 
/>Q2 = (1— 6^0/(1 is the probability of fixation of a 

migrant with genotype '2' in the large population, which is 
assumed to exhibit genotype '0' before migration (see Eq. 1). is 
the average of the minimum crossing time among S independent 
islands. We again focus, for simplicity, on the limit where the first 
step of valley or plateau crossing, which occurs at rate roi, is much 
longer than the second. Then, we simply have Zc^'^ii/'S (see 
Results). In this expression, zaKr^^^ = l/{Nndpoi) (with poi 
obtained from Eq. 1) is the average crossing time for an isolated 
island. Hence, the champion island crosses the valley S times faster 
on average than a single isolated island. The second condition, 

< Zc, finally yields 



S^Npoi M 



(48) 



Together, Eqs. 47 and 48 yield the interval of M/ fid over which 
we expect subdivision to maximally accelerate crossing: 



S^Npoi M NSpn 
— 1 — « « 



(49) 



In this range, we expect the valley or plateau crossing time T/ of 
the large population to be dominated by the crossing time of the 
champion island, so that t//t// % 1 /S. This prediction is confirmed 
by our simulations (see Fig. 4B). 
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